﻿/* RDCorrection.h
 * 
 * Copyright (C) 2020, 2021 Yuan LI
 * email:18069211@qq.com
 * https://blog.csdn.net/liyuanbhu
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef RDCORRECTION_H
#define RDCORRECTION_H

#include <QImage>

/**
 * @brief The RDCorrection class 镜头的畸变修正
 * 只考虑了镜头的径向畸变，也就是桶形畸变和枕形畸变。
 * 使用的模型比较简单，只有一个参数 k > 0 修正桶形畸变，k < 0 修正枕形畸变。
 * 修正后的图像会有一定的放缩。四周也有可能出现黑边或者切掉一些。
 * 这里假设径向畸变可以表示为：
 * \bar u = \frac{2}{1+\sqrt{1- 4k (u^2+v^2)}} u
 * \bar v = \frac{2}{1+\sqrt{1- 4k (u^2+v^2)}} v
 * 那么修正公式为：
 * u = \frac{1}{1 + k (\bar u^2 + \bar v^2)} \bar u
 * v = \frac{1}{1 + k (\bar u^2 + \bar v^2)} \bar v
 * 这个公式参考了 Carsten Steger 等所著 Machine Vision Algorithms and Applications 一书 3.9 节。
 */
class RDCorrection
{
public:
    RDCorrection();

    /**
     * @brief map_nearest 镜头畸变修正，输出图像与输入图像同样大小。采用最临近插值
     * @param in 输入的图像，灰度图或者 RGB32 类型
     * @param k 修正系数，k > 0 修正桶形畸变，k < 0 修正枕形畸变。
     * @return 修正后的图像
     */
    QImage map_nearest(const QImage &in, double k);
    /**
     * @brief map_nearest 镜头畸变修正，输出图像与输入图像同样大小。采用双线性插值
     * @param in 输入的图像，灰度图或者 RGB32 类型
     * @param k 修正系数，k > 0 修正桶形畸变，k < 0 修正枕形畸变。
     * @return 修正后的图像
     */
    QImage map_bilinear(const QImage &in, double k);
    /**
     * @brief map_nearest 镜头畸变修正，采用最临近插值。输出的图像可以与输入图像不相同的尺寸。
     * @param in 输入图像。
     * @param k 修正系数，k > 0 修正桶形畸变，k < 0 修正枕形畸变。
     * @param out 输出图像。传入这个参数时要求 format 与输入图像相同。
     * @return true 表示成功。
     */
    bool map_nearest(const QImage &in, double k, QImage &out);
    /**
     * @brief map_nearest 镜头畸变修正，采用双线性插值。输出的图像可以与输入图像不相同的尺寸。
     * @param in 输入图像。
     * @param k 修正系数，k > 0 修正桶形畸变，k < 0 修正枕形畸变。
     * @param out 输出图像。传入这个参数时要求 format 与输入图像相同。
     * @return true 表示成功。
     */
    bool map_bilinear(const QImage &in, double k, QImage &out);
    void setBackground(QRgb bg);
private:
    bool map_rgb32_nearest(const QImage &in, double k, QImage &out);
    bool map_rgb32_bilinear(const QImage &in, double k, QImage &out);
    bool map_gray_nearest(const QImage &in, double k, QImage &out);
    bool map_gray_bilinear(const QImage &in, double k, QImage &out);
    QImage map_gray_nearest(const QImage &in, double k);
    QImage map_gray_bilinear(const QImage &in, double k);
    QImage map_rgb32_nearest(const QImage &in, double k);
    QImage map_rgb32_bilinear(const QImage &in, double k);

    inline uchar bilinear(uchar in00, uchar in01, uchar in10, uchar in11, double x, double y)
    {
        return in00 + (in01 - in00) * x + (in10 - in00) * y + (in11 - in01 - in10 + in00) * x * y;
    }

    inline QRgb bilinear(QRgb in00, QRgb in01, QRgb in10, QRgb in11, double x, double y)
    {
        // f= b1 + b2 * x + b3 * y + b4 * x * y
        // b1 = in(0,0)
        // b2 = in(1,0) - in(0,0)
        // b3 = in(0,1) - in(0,0)
        // b4 = in(1,1) - in(1,0) - in(0,1) + in(0,0)
        int c00 = qRed(in00);
        int c01 = qRed(in01);
        int c10 = qRed(in10);
        int c11 = qRed(in11);

        int red = c00 + (c01 - c00) * x + (c10 - c00) * y + (c11 - c01 - c10 + c00) * x * y;
//        int red = (1 - x) * (1 - y) * c00 + (1 - x) * y * c01 + (1 - y) * x * c10 + x * y * c11;

        c00 = qGreen(in00);
        c01 = qGreen(in01);
        c10 = qGreen(in10);
        c11 = qGreen(in11);

        int green = c00 + (c01 - c00) * x + (c10 - c00) * y + (c11 - c01 - c10 + c00) * x * y;//
        //int green = (1 - x) * (1 - y) * c00 + (1 - x) * y * c01 + (1 - y) * x * c10 + x * y * c11;

        c00 = qBlue(in00);
        c01 = qBlue(in01);
        c10 = qBlue(in10);
        c11 = qBlue(in11);

        int blue =  c00 + (c01 - c00) * x + (c10 - c00) * y + (c11 - c01 - c10 + c00) * x * y;
        //int blue = (1 - x) * (1 - y) * c00 + (1 - x) * y * c01 + (1 - y) * x * c10 + x * y * c11;

        return qRgb(red, green, blue);
    }

    QRgb m_bg_rgb32;
    uchar m_bg_gray;
};

#endif // RDCORRECTION_H
